Load in data
#load in global data set
globe_ice.df <- read.csv(file = "./data/global_sea_ice_1979-2022.csv", header = TRUE) %>%
as_tibble() %>%
rename("global_ice" = "Ice",
"global_anomaly" = "Anomaly")
summary(globe_ice.df)
## Date Date2 Year Month
## Min. :197901 Length:518 Min. :1979 Min. : 1.00
## 1st Qu.:198910 Class :character 1st Qu.:1989 1st Qu.: 3.00
## Median :200008 Mode :character Median :2000 Median : 6.00
## Mean :200015 Mean :2000 Mean : 6.48
## 3rd Qu.:201105 3rd Qu.:2011 3rd Qu.: 9.00
## Max. :202202 Max. :2022 Max. :12.00
## global_ice global_anomaly
## Min. :-9999 Min. :-9999
## 1st Qu.: 20 1st Qu.: -1
## Median : 23 Median : 0
## Mean : -16 Mean : -39
## 3rd Qu.: 25 3rd Qu.: 0
## Max. : 28 Max. : 2
#load in northern hemisphere data set
nohem_ice.df <- read.csv(file = "./data/nohem_sea_ice_1979-2022.csv", header = TRUE) %>%
as_tibble() %>%
rename("nohem_ice" = "Value",
"nohem_anomaly" = "Anomaly")
summary(nohem_ice.df)
## Date Date2 Year Month
## Min. :197901 Length:519 Min. :1979 Min. : 1.00
## 1st Qu.:198910 Class :character 1st Qu.:1989 1st Qu.: 3.00
## Median :200008 Mode :character Median :2000 Median : 6.00
## Mean :200019 Mean :2000 Mean : 6.47
## 3rd Qu.:201106 3rd Qu.:2011 3rd Qu.: 9.00
## Max. :202203 Max. :2022 Max. :12.00
## nohem_ice nohem_anomaly
## Min. :-9999 Min. :-9999
## 1st Qu.: 9 1st Qu.: -1
## Median : 12 Median : 0
## Mean : -27 Mean : -39
## 3rd Qu.: 14 3rd Qu.: 0
## Max. : 16 Max. : 1
#load in southern hemisphere data set
sohem_ice.df <- read.csv(file = "./data/sohem_sea_ice_1979-2022.csv", header = TRUE) %>%
as_tibble() %>%
rename("sohem_ice" = "Value",
"sohem_anomaly" = "Anomaly")
summary(sohem_ice.df)
## Date Date2 Year Month
## Min. :197901 Length:519 Min. :1979 Min. : 1.00
## 1st Qu.:198910 Class :character 1st Qu.:1989 1st Qu.: 3.00
## Median :200008 Mode :character Median :2000 Median : 6.00
## Mean :200019 Mean :2000 Mean : 6.47
## 3rd Qu.:201106 3rd Qu.:2011 3rd Qu.: 9.00
## Max. :202203 Max. :2022 Max. :12.00
## sohem_ice sohem_anomaly
## Min. :-9999 Min. :-9999
## 1st Qu.: 6 1st Qu.: 0
## Median : 12 Median : 0
## Mean : -27 Mean : -39
## 3rd Qu.: 17 3rd Qu.: 0
## Max. : 20 Max. : 2
#load in hemispheric daily sea ice extent for duration calculations/plots
nohem_ice_daily.df <- read.csv(file = "./data/N_seaice_extent_daily_v3.0.csv", header = TRUE) %>%
as_tibble() %>%
rename("nohem_extent" = "Extent") %>%
mutate(Date = as.Date(Date)) %>%
filter(!Date < "1988-12-31" & Year != 2022) #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation, also remove 2022 since incomplete data
summary(nohem_ice_daily.df)
## Date Year Month Day
## Min. :1988-12-31 Min. :1988 Min. : 1.00 Min. : 1.0
## 1st Qu.:1997-04-01 1st Qu.:1997 1st Qu.: 4.00 1st Qu.: 8.0
## Median :2005-07-01 Median :2005 Median : 7.00 Median :16.0
## Mean :2005-07-01 Mean :2005 Mean : 6.52 Mean :15.7
## 3rd Qu.:2013-09-30 3rd Qu.:2013 3rd Qu.:10.00 3rd Qu.:23.0
## Max. :2021-12-31 Max. :2021 Max. :12.00 Max. :31.0
## nohem_extent
## Min. : 3.34
## 1st Qu.: 8.30
## Median :11.86
## Mean :11.11
## 3rd Qu.:14.09
## Max. :16.25
sohem_ice_daily.df <- read.csv(file = "./data/S_seaice_extent_daily_v3.0.csv", header = TRUE) %>%
as_tibble() %>%
rename("sohem_extent" = "Extent") %>%
mutate(Date = as.Date(Date)) %>%
filter(!Date < "1988-12-31") %>% #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation
filter(!Date < "1988-12-31" & Year != 2022) #remove up to last day of 1987, years prior sampled every other day, need last day for lag calculation, also remove 2022 since incomplete data
summary(sohem_ice_daily.df)
## Date Year Month Day
## Min. :1988-12-31 Min. :1988 Min. : 1.00 Min. : 1.0
## 1st Qu.:1997-04-01 1st Qu.:1997 1st Qu.: 4.00 1st Qu.: 8.0
## Median :2005-07-01 Median :2005 Median : 7.00 Median :16.0
## Mean :2005-07-01 Mean :2005 Mean : 6.52 Mean :15.7
## 3rd Qu.:2013-09-30 3rd Qu.:2013 3rd Qu.:10.00 3rd Qu.:23.0
## Max. :2021-12-31 Max. :2021 Max. :12.00 Max. :31.0
## sohem_extent
## Min. : 2.08
## 1st Qu.: 6.09
## Median :12.59
## Mean :11.68
## 3rd Qu.:17.24
## Max. :20.20
#load in sea level data
sea_level.df <- read.csv(file = "./data/epa-sea-level.csv", header = TRUE) %>%
as_tibble() %>%
rename("mean_sea_level" = "CSIRO.Adjusted.Sea.Level",
"lowerCI_sea_level" = "Lower.Error.Bound",
"upperCI_sea_level" = "Upper.Error.Bound")
summary(sea_level.df)
## Date Year mean_sea_level lowerCI_sea_level
## Length:35 Min. :1979 Min. :5.36 Min. :5.10
## Class :character 1st Qu.:1988 1st Qu.:6.16 1st Qu.:5.91
## Mode :character Median :1996 Median :6.67 Median :6.39
## Mean :1996 Mean :6.99 Mean :6.72
## 3rd Qu.:2004 3rd Qu.:7.75 3rd Qu.:7.47
## Max. :2013 Max. :9.33 Max. :8.99
## upperCI_sea_level
## Min. :5.63
## 1st Qu.:6.41
## Median :6.94
## Mean :7.26
## 3rd Qu.:8.03
## Max. :9.66
#load in temp data
global_temp.df <- read.csv(file = "./data/monthly_globaltemp.csv", header = TRUE) %>%
as_tibble() %>%
arrange(-row_number()) %>% #reverse row order so earliest years come first, matches sea ice df
mutate(Date = as.Date(Date) %>%
as.yearmon()) #change to month year format
#join the three dataframes
total_ice.df <- inner_join(globe_ice.df, nohem_ice.df,
by = c("Date" = "Date",
"Date2" = "Date2",
"Year" = "Year",
"Month" = "Month"))
total_ice.df <- inner_join(total_ice.df, sohem_ice.df,
by = c("Date" = "Date",
"Date2" = "Date2",
"Year" = "Year",
"Month" = "Month"))
#tidy up data
total_ice.df <- total_ice.df %>%
filter(global_ice != -9999) %>% #remove rows with -9999 (no value)
dplyr::select(!Date) %>%
rename(Date = Date2) %>%
mutate(Date = str_replace_all(string = Date,
pattern = "/",
replacement = "-"))
total_ice.df <- total_ice.df %>%
mutate(Date = mdy(Date))
#check how summary looks now
summary(total_ice.df)
## Date Year Month global_ice
## Min. :1979-01-01 Min. :1979 Min. : 1.00 Min. :16.3
## 1st Qu.:1989-11-23 1st Qu.:1989 1st Qu.: 3.00 1st Qu.:20.3
## Median :2000-08-16 Median :2000 Median : 6.00 Median :23.5
## Mean :2000-08-03 Mean :2000 Mean : 6.48 Mean :23.0
## 3rd Qu.:2011-05-08 3rd Qu.:2011 3rd Qu.: 9.00 3rd Qu.:25.3
## Max. :2022-02-01 Max. :2022 Max. :12.00 Max. :27.8
## global_anomaly nohem_ice nohem_anomaly sohem_ice
## Min. :-3.720 Min. : 3.57 Min. :-3.020 Min. : 2.16
## 1st Qu.:-0.725 1st Qu.: 8.56 1st Qu.:-0.843 1st Qu.: 5.82
## Median :-0.100 Median :12.10 Median :-0.185 Median :12.12
## Mean :-0.279 Mean :11.39 Mean :-0.274 Mean :11.57
## 3rd Qu.: 0.390 3rd Qu.:14.30 3rd Qu.: 0.340 3rd Qu.:17.02
## Max. : 1.650 Max. :16.34 Max. : 1.260 Max. :19.76
## sohem_anomaly
## Min. :-2.1300
## 1st Qu.:-0.3200
## Median :-0.0350
## Mean :-0.0059
## 3rd Qu.: 0.3300
## Max. : 1.8500
str(total_ice.df)
## tibble [516 × 9] (S3: tbl_df/tbl/data.frame)
## $ Date : Date[1:516], format: "1979-01-01" "1979-02-01" ...
## $ Year : int [1:516] 1979 1979 1979 1979 1979 1979 1979 1979 1979 1979 ...
## $ Month : int [1:516] 1 2 3 4 5 6 7 8 9 10 ...
## $ global_ice : num [1:516] 20.8 19.3 20.3 22.9 24.7 ...
## $ global_anomaly: num [1:516] 1.39 0.95 0.88 1.4 1.24 1.6 1.4 0.82 0.35 0.13 ...
## $ nohem_ice : num [1:516] 15.4 16.2 16.3 15.4 13.9 ...
## $ nohem_anomaly : num [1:516] 0.99 0.88 0.91 0.76 0.57 0.76 0.84 0.84 0.64 0.4 ...
## $ sohem_ice : num [1:516] 5.4 3.14 4 7.49 10.83 ...
## $ sohem_anomaly : num [1:516] 0.4 0.07 -0.03 0.64 0.67 0.84 0.56 -0.02 -0.29 -0.28 ...
total_ice_bar.df <- total_ice.df %>%
pivot_longer(cols = ends_with("ice"), #taking from wide to long
names_to = "ice_type", #category = type of event
values_to = "area") %>% #measurement = time till event
filter(ice_type != "global_hem") %>%
mutate(ice_type = as.factor(ice_type) %>%
fct_recode("Northern" = "nohem_ice",
"Southern" = "sohem_ice"))
#make ggplotly
sea.ice.gg <- total_ice_bar.df %>%
ggplot(aes(x = Date, y = area, fill = ice_type)) +
geom_bar(stat = "identity", position = "stack", size = 1.5) +
scale_y_continuous(limits = c(0,25)) +
scale_x_date( #The date equiv of scale_x_continuous
date_breaks = "5 year", #Breaks
date_minor_breaks = "1 year", #Little lines in between
date_labels = "%b '%y" #Label for the date,
) +
xlab("") + ylab("Ice Area (millions km^2)") +
scale_fill_manual(values = met.brewer("Hokusai3", direction = 1))
#ggplotly() + tooltip to control label
ggplotly(sea.ice.gg, tooltip = "text") #tooltip controls text hover/label
total_ice_line.df <- total_ice.df %>%
pivot_longer(cols = ends_with("ice"), #taking from wide to long
names_to = "ice_type", #category = type of event
values_to = "extent") %>% #measurement = time till event
mutate(ice_type = as.factor(ice_type) %>%
fct_recode("Global" = "global_ice",
"Northern" = "nohem_ice",
"Southern" = "sohem_ice"))
Now make in plotly
#make plotly of sea ice
seaice.plotly <- plot_ly(total_ice_line.df, x = ~Date, y = ~extent, color = ~ice_type)
seaice.plotly <- seaice.plotly %>% add_lines()
seaice.plotly <- seaice.plotly %>% layout(
xaxis = list(
rangeselector = list(
buttons = list(
list(
count = 6,
label = "6 mo",
step = "month",
stepmode = "backward"),
list(
count = 1,
label = "1 yr",
step = "year",
stepmode = "backward"),
list(
count = 5,
label = "5 yr",
step = "year",
stepmode = "backward"),
list(
count = 10,
label = "10 yr",
step = "year",
stepmode = "backward"),
list(
count = 1,
label = "YTD",
step = "year",
stepmode = "todate"),
list(step = "all"))),
rangeslider = list(type = "date")),
yaxis = list(title = "Sea Ice Extent (millions km^2)"))
seaice.plotly
Let’s make some plots looking at annual, as opposed to monthly, shifts in sea ice
total_ice_year.df <- total_ice_line.df %>%
group_by(Year, ice_type) %>%
summarize(
mean_ice_extent = mean(extent)
) %>%
ungroup() %>%
filter(Year != 2022) #remove last year since incomplete, throws off annual mean
#now plot
seaice.year.plotly <- plot_ly(data = total_ice_year.df, x = ~Year, y = ~mean_ice_extent, color = ~ice_type)
seaice.year.plotly <- seaice.year.plotly %>% add_lines()
seaice.year.plotly <- seaice.year.plotly %>% layout(
xaxis = list("Year"),
yaxis = list(title = "Sea Ice Extent (millions km^2)"))
seaice.year.plotly
Not surprisingly, it looks like we have a declining trend in sea ice extent. Let’s create some simple linear models to compare trends.
#make global dataset
global_ice_year.df <- total_ice_year.df %>%
filter(ice_type %in% "Global")
#make northern dataset
nohem_ice_year.df <- total_ice_year.df %>%
filter(ice_type %in% "Northern")
#make southern dataset
sohem_ice_year.df <- total_ice_year.df %>%
filter(ice_type %in% "Southern")
#construct three simple linear models comparing year with ice extent to get annual mean trend
mod1 <- lm(mean_ice_extent ~ Year, data = global_ice_year.df)
mod2 <- lm(mean_ice_extent ~ Year, data = nohem_ice_year.df)
mod3 <- lm(mean_ice_extent ~ Year, data = sohem_ice_year.df)
#check out model summaries
summary(mod1)
##
## Call:
## lm(formula = mean_ice_extent ~ Year, data = global_ice_year.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.0262 -0.3830 -0.0058 0.3960 1.2481
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 121.21108 12.45421 9.73 0.0000000000032 ***
## Year -0.04912 0.00623 -7.89 0.0000000009479 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.507 on 41 degrees of freedom
## Multiple R-squared: 0.603, Adjusted R-squared: 0.593
## F-statistic: 62.2 on 1 and 41 DF, p-value: 0.000000000948
summary(mod2)
##
## Call:
## lm(formula = mean_ice_extent ~ Year, data = nohem_ice_year.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5010 -0.1687 0.0225 0.1707 0.3300
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 119.55461 5.33322 22.4 <0.0000000000000002 ***
## Year -0.05409 0.00267 -20.3 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.217 on 41 degrees of freedom
## Multiple R-squared: 0.909, Adjusted R-squared: 0.907
## F-statistic: 411 on 1 and 41 DF, p-value: <0.0000000000000002
summary(mod3)
##
## Call:
## lm(formula = mean_ice_extent ~ Year, data = sohem_ice_year.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.9825 -0.2006 0.0013 0.1675 1.0565
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.67151 10.20383 0.16 0.87
## Year 0.00496 0.00510 0.97 0.34
##
## Residual standard error: 0.415 on 41 degrees of freedom
## Multiple R-squared: 0.0226, Adjusted R-squared: -0.00127
## F-statistic: 0.947 on 1 and 41 DF, p-value: 0.336
Based on these simple linear models, if mean annual sea ice trends continued at the same rate as they have since 1979, we would predict that average sea ice extent will reach 0 in the northern hemisphere by the year 2210, and it will reach 0 globally and, thus, in the southern hemisphere, by the year 2468.
Let’s calculate the duration of sea ice expansion and and recession for each year
#join both daily data frames
bothhem_ice_daily.df <- inner_join(nohem_ice_daily.df, sohem_ice_daily.df,
by = c("Date" = "Date",
"Year" = "Year",
"Month" = "Month",
"Day" = "Day"))
#for northern hemisphere
bothhem_ice_daily.df <- bothhem_ice_daily.df %>%
mutate(nohem_ice_diff = nohem_extent - lag(nohem_extent),
sohem_ice_diff = sohem_extent - lag(sohem_extent),
nohem_ice_trend = case_when(
nohem_ice_diff > 0 ~ "increasing",
nohem_ice_diff %in% 0 ~ "no change",
nohem_ice_diff < 0 ~ "decreasing"),
sohem_ice_trend = case_when(
sohem_ice_diff > 0 ~ "increasing",
sohem_ice_diff %in% 0 ~ "no change",
sohem_ice_diff < 0 ~ "decreasing")
) %>%
na.omit()
summary(bothhem_ice_daily.df)
## Date Year Month Day
## Min. :1989-01-01 Min. :1989 Min. : 1.00 Min. : 1.0
## 1st Qu.:1997-04-02 1st Qu.:1997 1st Qu.: 4.00 1st Qu.: 8.0
## Median :2005-07-02 Median :2005 Median : 7.00 Median :16.0
## Mean :2005-07-02 Mean :2005 Mean : 6.52 Mean :15.7
## 3rd Qu.:2013-10-01 3rd Qu.:2013 3rd Qu.:10.00 3rd Qu.:23.0
## Max. :2021-12-31 Max. :2021 Max. :12.00 Max. :31.0
## nohem_extent sohem_extent nohem_ice_diff sohem_ice_diff
## Min. : 3.34 Min. : 2.08 Min. :-0.3970 Min. :-0.5180
## 1st Qu.: 8.30 1st Qu.: 6.09 1st Qu.:-0.0570 1st Qu.:-0.0730
## Median :11.86 Median :12.59 Median :-0.0050 Median : 0.0200
## Mean :11.11 Mean :11.68 Mean :-0.0001 Mean :-0.0001
## 3rd Qu.:14.09 3rd Qu.:17.24 3rd Qu.: 0.0540 3rd Qu.: 0.0880
## Max. :16.25 Max. :20.20 Max. : 0.3630 Max. : 0.3310
## nohem_ice_trend sohem_ice_trend
## Length:12053 Length:12053
## Class :character Class :character
## Mode :character Mode :character
##
##
##
#now summarize number of increasing and decreasing days for each year
#northern hemisphere
nohem_ice_duration.df <- bothhem_ice_daily.df %>%
group_by(Year, nohem_ice_trend) %>%
summarise(
no_days = length(nohem_ice_trend)
) %>%
mutate(nohem_ice_trend = as.factor(nohem_ice_trend)) %>%
ungroup()
#southern hemisphere
sohem_ice_duration.df <- bothhem_ice_daily.df %>%
group_by(Year, sohem_ice_trend) %>%
summarise(
no_days = length(sohem_ice_trend)
) %>%
mutate(sohem_ice_trend = as.factor(sohem_ice_trend)) %>%
ungroup()
Now make plots showing trends in days if increasing and decreasing ice
#northern hemisphere plot
nohem.ice.dur.plotly <- plot_ly(data = nohem_ice_duration.df, x = ~Year, y = ~no_days, color = ~nohem_ice_trend)
nohem.ice.dur.plotly <- nohem.ice.dur.plotly %>% add_lines()
nohem.ice.dur.plotly <- nohem.ice.dur.plotly %>% layout(
xaxis = list("Year"),
yaxis = list(title = "No. of Days"))
nohem.ice.dur.plotly
#southern hemisphere plot
sohem.ice.dur.plotly <- plot_ly(data = sohem_ice_duration.df, x = ~Year, y = ~no_days, color = ~sohem_ice_trend)
sohem.ice.dur.plotly <- sohem.ice.dur.plotly %>% add_lines()
sohem.ice.dur.plotly <- sohem.ice.dur.plotly %>% layout(
xaxis = list("Year"),
yaxis = list(title = "No. of Days"))
sohem.ice.dur.plotly
Now do the same but for anomaly data
#alter df for anomaly plot
total_ice_anomaly.df <- total_ice.df %>%
pivot_longer(cols = ends_with("anomaly"), #taking from wide to long
names_to = "ice_type", #category = type of event
values_to = "anomaly") %>% #measurement = time till event
mutate(ice_type = as.factor(ice_type) %>%
fct_recode("Global" = "global_anomaly",
"Northern" = "nohem_anomaly",
"Southern" = "sohem_anomaly"))
#make plotly of sea ice anomaly values
iceanomaly.plotly <- plot_ly(total_ice_anomaly.df, x = ~Date, y = ~anomaly, color = ~ice_type)
iceanomaly.plotly <- iceanomaly.plotly %>% add_lines()
iceanomaly.plotly <- iceanomaly.plotly %>% layout(
xaxis = list(
rangeselector = list(
buttons = list(
list(
count = 6,
label = "6 mo",
step = "month",
stepmode = "backward"),
list(
count = 1,
label = "1 yr",
step = "year",
stepmode = "backward"),
list(
count = 5,
label = "5 yr",
step = "year",
stepmode = "backward"),
list(
count = 10,
label = "10 yr",
step = "year",
stepmode = "backward"),
list(
count = 1,
label = "YTD",
step = "year",
stepmode = "todate"),
list(step = "all"))),
rangeslider = list(type = "date")),
yaxis = list(title = "Sea Ice Anomalies"))
iceanomaly.plotly
It looks like global anomalies in sea ice are mostly driven by the northern hemisphere sea ice as opposed to the southern hemisphere. Let’s rtun some simple linear models to check it out.
#global anomaly vs southern hemisphere
mod4 <- lm(global_anomaly ~ sohem_anomaly, data = total_ice.df)
summary(mod4)
##
## Call:
## lm(formula = global_anomaly ~ sohem_anomaly, data = total_ice.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.692 -0.546 0.113 0.594 1.571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.2739 0.0348 -7.87 0.000000000000021 ***
## sohem_anomaly 0.8863 0.0625 14.19 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.791 on 514 degrees of freedom
## Multiple R-squared: 0.281, Adjusted R-squared: 0.28
## F-statistic: 201 on 1 and 514 DF, p-value: <0.0000000000000002
#global anomaly vs northern hemisphere
mod5 <- lm(global_anomaly ~ nohem_anomaly, data = total_ice.df)
summary(mod5)
##
## Call:
## lm(formula = global_anomaly ~ nohem_anomaly, data = total_ice.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.1877 -0.3161 -0.0229 0.3070 1.8240
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0205 0.0259 -0.79 0.43
## nohem_anomaly 0.9434 0.0309 30.49 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.557 on 514 degrees of freedom
## Multiple R-squared: 0.644, Adjusted R-squared: 0.643
## F-statistic: 930 on 1 and 514 DF, p-value: <0.0000000000000002
#plots of both models above
#sohem
total_ice.df %>%
ggplot(aes(x = sohem_anomaly, y = global_anomaly)) +
stat_smooth(formula = "y ~ x",
method = "lm")
#nohem
total_ice.df %>%
ggplot(aes(x = nohem_anomaly, y = global_anomaly)) +
stat_smooth(formula = "y ~ x",
method = "lm")
As expected from an earlier visual assessment of anomaly data, the model with northern hemisphere anomaly explains 64.3% of variation in global anomaly data, relative to only 28% explained by the southern hemisphere anomalies.
#prep sea ice df for comparison with sea level data
total_ice_march.df <- total_ice.df %>%
filter(Month == 3) %>% #only use ice data from March, matches collectiion of annual sea level data
dplyr::select(Year, global_ice, global_anomaly) %>%
mutate(global_ice_change = global_ice - lag(global_ice)) %>% #create annual change variable
filter(Year != 1979) #remove first year with missing change value
#remove Date column
sea_level.df <- sea_level.df %>%
dplyr::select(!Date)
#join dataframes
sea_levelandice.df <- inner_join(total_ice_march.df, sea_level.df,
by = ("Year" = "Year"))
#check out summary
summary(sea_levelandice.df)
## Year global_ice global_anomaly global_ice_change
## Min. :1980 Min. :17.6 Min. :-1.8300 Min. :-1.1500
## 1st Qu.:1988 1st Qu.:19.1 1st Qu.:-0.3375 1st Qu.:-0.5875
## Median :1996 Median :19.4 Median :-0.0100 Median :-0.1950
## Mean :1996 Mean :19.4 Mean :-0.0253 Mean :-0.0085
## 3rd Qu.:2005 3rd Qu.:19.9 3rd Qu.: 0.4525 3rd Qu.: 0.4750
## Max. :2013 Max. :20.5 Max. : 1.0000 Max. : 2.0900
## mean_sea_level lowerCI_sea_level upperCI_sea_level
## Min. :5.60 Min. :5.34 Min. :5.85
## 1st Qu.:6.17 1st Qu.:5.92 1st Qu.:6.42
## Median :6.73 Median :6.46 Median :7.00
## Mean :7.04 Mean :6.76 Mean :7.30
## 3rd Qu.:7.75 3rd Qu.:7.48 3rd Qu.:8.04
## Max. :9.33 Max. :8.99 Max. :9.66
#manual colors for legend
colors_icelevel <- c("Global Sea Level" = "#3399FF", "Global Sea Ice" = "#993399")
#plot sea ice and level together
sea_levelandice.df %>%
ggplot(aes(x = Year)) +
geom_line(aes(y = global_ice, color = "Global Sea Ice")) +
geom_line(aes(y = mean_sea_level*3, color = "Global Sea Level")) +
geom_ribbon(aes(ymin = lowerCI_sea_level*3, ymax = upperCI_sea_level*3), fill = "#99CCFF", alpha = 0.25) +
#add themes to clean up graph
scale_y_continuous(
# Features of the first axis
name = "Global Sea Ice Extent (millions km^2)",
# Add a second axis and specify its features
sec.axis = sec_axis(trans = ~./3, name = "Mean Global Sea Level Change (in)")) +
scale_color_manual("", #manual legend
values = colors_icelevel) + #set of colors matched with plot type
theme(legend.position = "bottom",
legend.key = element_rect(fill = "white"))
#construct simple linear regression models
mod6 <- lm(mean_sea_level ~ global_ice, data = sea_levelandice.df)
mod7 <- lm(mean_sea_level ~ global_ice_change + poly(global_ice_change, 2), data = sea_levelandice.df)
#model summary
summary(mod6)
##
## Call:
## lm(formula = mean_sea_level ~ global_ice, data = sea_levelandice.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.484 -0.823 -0.221 0.589 2.437
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 15.868 5.681 2.79 0.0087 **
## global_ice -0.454 0.292 -1.56 0.1297
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.07 on 32 degrees of freedom
## Multiple R-squared: 0.0703, Adjusted R-squared: 0.0412
## F-statistic: 2.42 on 1 and 32 DF, p-value: 0.13
summary(mod7)
##
## Call:
## lm(formula = mean_sea_level ~ global_ice_change + poly(global_ice_change,
## 2), data = sea_levelandice.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.859 -0.656 -0.208 0.559 2.352
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.037 0.170 41.31 <0.0000000000000002 ***
## global_ice_change 0.116 0.210 0.55 0.5839
## poly(global_ice_change, 2)1 NA NA NA NA
## poly(global_ice_change, 2)2 2.866 0.993 2.89 0.0071 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.993 on 31 degrees of freedom
## Multiple R-squared: 0.218, Adjusted R-squared: 0.167
## F-statistic: 4.32 on 2 and 31 DF, p-value: 0.0222
#check residuals vs fitted
par(mfrow = c(2,2))
plot(mod6)
par(mfrow = c(2,2))
plot(mod7)
#plot relationships
sea_levelandice.df %>%
ggplot(aes(x = global_ice, y = mean_sea_level)) +
stat_smooth(formula = "y ~ poly(x,2)",
method = "lm")
sea_levelandice.df %>%
ggplot(aes(x = global_ice_change, y = mean_sea_level)) +
stat_smooth(formula = "y ~ poly(x,2)",
method = "lm")
#convert date to month year to match temperature data
total_ice_yearmon.df <- total_ice.df %>%
mutate(Date = as.yearmon(Date))
#split temp df into two dfsbased on source
global_temp_gis.df <- global_temp.df %>%
filter(Source %in% "GISTEMP")
global_temp_gcag.df <- global_temp.df %>%
filter(Source %in% "GCAG")
#join dataframes
ice_temp_gis.df <- inner_join(total_ice_yearmon.df, global_temp_gis.df,
by = ("Date" = "Date"))
ice_temp_gcag.df <- inner_join(total_ice_yearmon.df, global_temp_gcag.df,
by = ("Date" = "Date"))
#manual colors for legend
colors_icetemp <- c("Global Temperature" = "#3399FF", "Global Sea Ice" = "#993399")
#plot sea ice and level together
ice_temp_gis.df %>%
ggplot(aes(x = Date)) +
geom_line(aes(y = global_ice, color = "Global Sea Ice")) +
geom_line(aes(y = Mean*25, color = "Global Temperature")) +
#add themes to clean up graph
scale_y_continuous(
# Features of the first axis
name = "Global Sea Ice Extent (millions km^2)",
# Add a second axis and specify its features
sec.axis = sec_axis(trans = ~./25, name = "Mean Global Temp Change (degrees Celsius)")) +
scale_color_manual("", #manual legend
values = colors_icetemp) + #set of colors matched with plot type
theme(legend.position = "bottom",
legend.key = element_rect(fill = "white"))
#plot sea ice and level together
ice_temp_gis.df %>%
ggplot(aes(x = Date)) +
geom_line(aes(y = global_anomaly, color = "Global Sea Ice")) +
geom_line(aes(y = Mean*3, color = "Global Temperature")) +
#add themes to clean up graph
scale_y_continuous(
# Features of the first axis
name = "Global Sea Ice Anomalies",
# Add a second axis and specify its features
sec.axis = sec_axis(trans = ~./3, name = "Mean Global Temp Change (degrees Celsius)")) +
scale_color_manual("", #manual legend
values = colors_icetemp) + #set of colors matched with plot type
theme(legend.position = "bottom",
legend.key = element_rect(fill = "white"))
#construct simple linear regression models
mod8 <- lm(global_ice ~ Mean, data = ice_temp_gis.df)
mod9 <- lm(global_anomaly ~ Mean, data = ice_temp_gis.df)
#model summary
summary(mod8)
##
## Call:
## lm(formula = global_ice ~ Mean, data = ice_temp_gis.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.227 -2.568 0.848 2.146 4.510
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 24.383 0.287 85.01 < 0.0000000000000002 ***
## Mean -2.571 0.546 -4.71 0.0000033 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.75 on 452 degrees of freedom
## Multiple R-squared: 0.0468, Adjusted R-squared: 0.0447
## F-statistic: 22.2 on 1 and 452 DF, p-value: 0.00000327
summary(mod9)
##
## Call:
## lm(formula = global_anomaly ~ Mean, data = ice_temp_gis.df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.8921 -0.3729 0.0605 0.4468 1.6615
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.6942 0.0694 10.0 <0.0000000000000002 ***
## Mean -1.6366 0.1319 -12.4 <0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.666 on 452 degrees of freedom
## Multiple R-squared: 0.254, Adjusted R-squared: 0.252
## F-statistic: 154 on 1 and 452 DF, p-value: <0.0000000000000002
#check residuals vs fitted
par(mfrow = c(2,2))
plot(mod8)
par(mfrow = c(2,2))
plot(mod9)
#plot relationships
ice_temp_gis.df %>%
ggplot(aes(x = Mean, y = global_ice)) +
stat_smooth(formula = "y ~ x",
method = "lm")
ice_temp_gis.df %>%
ggplot(aes(x = Mean, y = global_anomaly)) +
stat_smooth(formula = "y ~ x",
method = "lm")
Read in a few csvs to get started
#load shapefiles
seaice.jan1980 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp")
## Reading layer `extent_N_198001_polygon_v3.0' from data source
## `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 143 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -3400000 ymin: -4875000 xmax: 3600000 ymax: 5850000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
seaice.jan2020 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202001_polygon_v3.0.shp")
## Reading layer `extent_N_202001_polygon_v3.0' from data source
## `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202001_polygon_v3.0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 158 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -3375000 ymin: -4925000 xmax: 3600000 ymax: 5850000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
seaice.sep1980 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198009_polygon_v3.0.shp")
## Reading layer `extent_N_198009_polygon_v3.0' from data source
## `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198009_polygon_v3.0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 118 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -2550000 ymin: -3250000 xmax: 2350000 ymax: 2100000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
seaice.sep2020 <- st_read("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202009_polygon_v3.0.shp")
## Reading layer `extent_N_202009_polygon_v3.0' from data source
## `/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_202009_polygon_v3.0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 83 features and 1 field
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: -2525000 ymin: -2800000 xmax: 2350000 ymax: 4950000
## Projected CRS: NSIDC Sea Ice Polar Stereographic North
# seaice.jan1980.transform <- stTransform(seaice.jan1980, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
#make basic plot of layers
ggplot() +
geom_sf(data = seaice.jan1980, fill = "blue") +
geom_sf(data = seaice.jan2020, fill = "purple")
ggplot() +
geom_sf(data = seaice.sep1980, fill = "red") +
geom_sf(data = seaice.sep2020, fill = "green")
ggplot() +
geom_sf(data = seaice.jan2020, fill = "brown") +
geom_sf(data = seaice.sep2020, fill = "yellow")
#load shapefile
seaice.jan1980 <- readOGR("./data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp")
## OGR data source with driver: ESRI Shapefile
## Source: "/Users/jackrabe/Documents/PUBH 7462 Advanced Programming in R/Homework/mn-ice-trends/data/Monthly Sea Ice Extent/all_spatial_files/extent_N_198001_polygon_v3.0.shp", layer: "extent_N_198001_polygon_v3.0"
## with 143 features
## It has 1 fields
## Integer64 fields read as strings: FID
#check class and crs
class(seaice.jan1980)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
crs(seaice.jan1980)
## Coordinate Reference System:
## Deprecated Proj.4 representation:
## +proj=stere +lat_0=90 +lat_ts=70 +lon_0=-45 +x_0=0 +y_0=0 +a=6378273
## +rf=298.279411123064 +units=m +no_defs +type=crs
## WKT2 2019 representation:
## PROJCRS["NSIDC Sea Ice Polar Stereographic North",
## BASEGEOGCRS["Unspecified datum based upon the Hughes 1980 ellipsoid",
## DATUM["Not specified (based on Hughes 1980 ellipsoid)",
## ELLIPSOID["Hughes 1980",6378273,298.279411123064,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4054]],
## CONVERSION["US NSIDC Sea Ice polar stereographic north",
## METHOD["Polar Stereographic (variant B)",
## ID["EPSG",9829]],
## PARAMETER["Latitude of standard parallel",70,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8832]],
## PARAMETER["Longitude of origin",-45,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8833]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting (X)",south,
## MERIDIAN[45,
## ANGLEUNIT["degree",0.0174532925199433]],
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",south,
## MERIDIAN[135,
## ANGLEUNIT["degree",0.0174532925199433]],
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["Polar research."],
## AREA["Northern hemisphere - north of 60°N onshore and offshore, including Arctic."],
## BBOX[60,-180,90,180]],
## ID["EPSG",3411]]
extent(seaice.jan1980)
## class : Extent
## xmin : -3400000
## xmax : 3600000
## ymin : -4875000
## ymax : 5850000
seaice.jan1980.WGS84 <- spTransform(seaice.jan1980, CRS("+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"))
crs(seaice.jan1980.WGS84)
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs
## WKT2 2019 representation:
## GEOGCRS["unknown",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]],
## ID["EPSG",6326]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8901]],
## CS[ellipsoidal,2],
## AXIS["longitude",east,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]],
## AXIS["latitude",north,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433,
## ID["EPSG",9122]]]]
extent(seaice.jan1980.WGS84)
## class : Extent
## xmin : -180
## xmax : 179.7
## ymin : 38.01
## ymax : 83.67
broom::tidy(seaice.jan1980.WGS84)
## Regions defined for each Polygons
## # A tibble: 4,470 × 7
## long lat order hole piece group id
## <dbl> <dbl> <int> <lgl> <fct> <fct> <chr>
## 1 125. 38.8 1 FALSE 1 0.1 0
## 2 125. 38.7 2 FALSE 1 0.1 0
## 3 124. 39.4 3 FALSE 1 0.1 0
## 4 125. 39.5 4 FALSE 1 0.1 0
## 5 125. 39.5 5 FALSE 1 0.1 0
## 6 125. 39.7 6 FALSE 1 0.1 0
## 7 125. 39.7 7 FALSE 1 0.1 0
## 8 125. 39.5 8 FALSE 1 0.1 0
## 9 125. 39.6 9 FALSE 1 0.1 0
## 10 125. 39.2 10 FALSE 1 0.1 0
## # … with 4,460 more rows
Make basemap
#manual bounding box
seaice.box <- matrix(c(-180, 180, 38.01, 83.67),
ncol = 2,
byrow = TRUE)
colnames(seaice.box) <- c("min", "max")
rownames(seaice.box) <- c("x", "y")
seaice.box
## min max
## x -180.00 180.00
## y 38.01 83.67
#Get the base map (foundational layer)
seaice_base.map <- get_map(
location = seaice.box,
source = "stamen",
maptype = "watercolor",
crop = TRUE
)
## Source : http://tile.stamen.com/terrain/3/0/0.png
## Source : http://tile.stamen.com/terrain/3/1/0.png
## Source : http://tile.stamen.com/terrain/3/2/0.png
## Source : http://tile.stamen.com/terrain/3/3/0.png
## Source : http://tile.stamen.com/terrain/3/4/0.png
## Source : http://tile.stamen.com/terrain/3/5/0.png
## Source : http://tile.stamen.com/terrain/3/6/0.png
## Source : http://tile.stamen.com/terrain/3/7/0.png
## Source : http://tile.stamen.com/terrain/3/8/0.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/0.png.
## Source : http://tile.stamen.com/terrain/3/0/1.png
## Source : http://tile.stamen.com/terrain/3/1/1.png
## Source : http://tile.stamen.com/terrain/3/2/1.png
## Source : http://tile.stamen.com/terrain/3/3/1.png
## Source : http://tile.stamen.com/terrain/3/4/1.png
## Source : http://tile.stamen.com/terrain/3/5/1.png
## Source : http://tile.stamen.com/terrain/3/6/1.png
## Source : http://tile.stamen.com/terrain/3/7/1.png
## Source : http://tile.stamen.com/terrain/3/8/1.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/1.png.
## Source : http://tile.stamen.com/terrain/3/0/2.png
## Source : http://tile.stamen.com/terrain/3/1/2.png
## Source : http://tile.stamen.com/terrain/3/2/2.png
## Source : http://tile.stamen.com/terrain/3/3/2.png
## Source : http://tile.stamen.com/terrain/3/4/2.png
## Source : http://tile.stamen.com/terrain/3/5/2.png
## Source : http://tile.stamen.com/terrain/3/6/2.png
## Source : http://tile.stamen.com/terrain/3/7/2.png
## Source : http://tile.stamen.com/terrain/3/8/2.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/2.png.
## Source : http://tile.stamen.com/terrain/3/0/3.png
## Source : http://tile.stamen.com/terrain/3/1/3.png
## Source : http://tile.stamen.com/terrain/3/2/3.png
## Source : http://tile.stamen.com/terrain/3/3/3.png
## Source : http://tile.stamen.com/terrain/3/4/3.png
## Source : http://tile.stamen.com/terrain/3/5/3.png
## Source : http://tile.stamen.com/terrain/3/6/3.png
## Source : http://tile.stamen.com/terrain/3/7/3.png
## Source : http://tile.stamen.com/terrain/3/8/3.png
## Not Found (HTTP 404). Failed to aquire tile /terrain/3/8/3.png.
ggmap(seaice_base.map) +
geom_polygon(data = seaice.jan1980.WGS84, aes(x = long, y = lat))
## Regions defined for each Polygons
leaflet(data = seaice.jan1980.WGS84) %>%
addProviderTiles("CartoDB.Positron") %>%
addPolygons()